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We develop a numerical approach for quantifying entanglement in mixed quantum states by 
convex-roof entanglement measures, based on the optimal entanglement witness operator and the 
minimax optimization method. Our new approach is applicable to general entanglement measures 
and states, and it is an efficient alternative to the conventional approach based upon the optimal 
pure-state decomposition. Compared with the conventional one, it has two important merits (i) 
that the global optimality of the solution is quantitatively verifiable, and (ii) that the optimization 
is considerably simplified by exploiting the common symmetry of the target state and measure. To 
demonstrate the merits, we quantify Greenberger-Horne-Zeilinger (GHZ) entanglement in a class of 
three-qubit full-rank mixed states composed of the GHZ state, the W state, and the white noise, the 
simplest mixtures of states with different genuine multipartite entanglement, which have not been 
quantified before this work. We discuss some general properties of the form of the optimal witness 
operator, and of the convex structure of mixed states, which are related to the symmetry and the 
rank of states. 
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I. INTRODUCTION 

Quantum entanglement, a quantum correlation among 
subsystems, is at the core of strange quantum phenom- 
ena, and it is the key ingredient of quantum information 
processing pQ. Its detection and quantification has at- 
tracted much attention of varieties of scientific commu- 
nities for both fundamental and practical purposes. 

Entanglement measure £ (p) quantifies the entangle- 
ment contained in a quantum state p [2J. The defini- 
tion of £(p) is not unique and depends on the type of 
entanglement to be quantified, as there exist many dif- 
ferent types of entanglement. A measure is constructed 
under the restrictions by the physical meaning of entan- 
glement: (i) Any measure vanishes for separable states 
p = J2iPi PA,i ® PB,i (8 • • • j (ii) is invariant under lo- 
cal unitary operations, and (iii) does not increase under 
local quantum operations and classical communications 
(LOCC). The last property, so-called monotonicity under 
LOCC, is usually treated by adopting the convexity of a 
measure, £(pp 1 + (l-p)p 2 ) < p£ (pi) + (l-p)£ (p 2 ) 0- A 
convex-roof measure, the most popular type of measure, 
is defined first for pure states \ip), and then extended, 
based on the convexity, to mixed states p via the convex 
roof construction 

£(p)= inf 5>£(hfc», (1) 
{pi,if>i} . 

where the optimization ("inf") runs over every possible 
pure state decomposition of p = Y^iPi n i>i an d = 
denotes the pure-state projector. 
In most cases, except for two-qubit states [1HS], 
the quantification of mixed-state entanglement requires 
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heavy computational cost for the optimization of a large 
number of parameters in Eq. [S] . The pure-state 
decomposition of a rank-n mixed state into m pure states 
utilizes a mxn left-unitary matrix which has 2mn— n 2 — 1 
parameters. As m ranges from n to n 2 according to 
Caratheodory's theorem [5], the maximal number of pa- 
rameters for the optimization becomes ~ 2n 3 ; it already 
reaches a thousand for three-qubit full-rank states. The 
optimization of such a large number of parameters is 
prone to the convergence to a local optimum, resulting in 
the overestimation of £(p). In usual optimization prob- 
lems, it cannot be ensured, without the overall investi- 
gation of the parameter space, whether the result of the 
optimization is globally optimal, i.e. provides the right 
answer. These facts make the quantification of multipar- 
tite and/or high-dimensional mixed-state entanglement 
practically unaccessible theoretically and experimentally 

There exists another definition of £{p), which is equiv- 
alent to Eq. 0. It is the dual problem of Eq. 0) in 
the sense that it is based on the optimization of quan- 
tum operators rather than the optimization of states; see 
Eq. ([2]) below. The optimized operator, whose expecta- 
tion value provides £(p), is called the optimal witness 
operator [TUlI12| ; here the optimal witness operator is 
used for entanglement quantification rather than detec- 
tion P21 [14] . This concept has been proved to be useful 
for estimating a lower bound of £{p) with partial infor- 
mation of p accessed in experiments [15j . 

Very recently, it was proposed |16j by two of the au- 
thors of the present work that an efficient way of con- 
structing the optimal witness operator can relieve the 
difficulty of optimizing £{p). The benefits of this ap- 
proach as an optimization problem are: (i) The size of 
the parameter space for this optimization is much smaller 
than that for Eq. , and (ii) one can verify whether the 
result of this approach is globally optimal, i.e. provides 
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the true answer, based on the duality between this ap- 
proach and Eq. (!]). This approach is generally applica- 
ble to any £ and p. This approach successfully quanti- 
fies thrcc-qubit Grccnbcrgcr-Hornc-Zcilingcr (GHZ) en- 
tanglement and four-qubit bound entanglement in some 
full-rank mixed states p that have been usually prepared 
in laboratories, even providing the analytic expression 
of £ (p), demonstrating its usefulness [16] ; this was, to 
our knowledge, the first quantification of entanglement 
in multipartite full-rank mixed states. 

In this paper, we formulate how to numerically obtain 
the optimal witness operator, utilizing the minimax opti- 
mization method [17] ; the utilization is shortly mentioned 
in Rcf. [T^] , but its detail was not investigated. To apply 
the minimax method, we derive the canonical form of the 
witness operator that is generally applicable to any mea- 
sures and states; see Theorem 1 below. We devise the 
techniques for handling the non-analyticity in the min- 
imax method and of reinforcing the robustness against 
numerical errors, and also exploit the existing theories of 
the minimax method [171 118] to enhance computation ef- 
ficiency. We emphasize that in our approach, the global 
optimality of a witness operator is quantitatively verified 
by the minimal distance between the target state and the 
set of the states exactly quantified by the witness opera- 
tor. We also notice that the common symmetries of the 
target state and measure simplify the optimization. 

Our method is applied to scrutinize three-qubit GHZ 
entanglement in exemplar states. In particular, we quan- 
tify GHZ entanglement in a family of full-rank mixed 
states composed of the GHZ state, the W state, and the 
white noise, which is the simplest full-rank mixture of 
pure states with different genuine multipartite entangle- 
ment; this type of mixed states has not been quantified 
before, to our knowledge. This reveals some general prop- 
erties of the form of the optimal witness operator, and of 
the convex structure of mixed states, which are related 
to the symmetry and the rank of the states. 

This paper is organized as follows. In S ec.|TT} we intro- 
duce the optimal witness operator. In Sec. |III[ we develop 
the minimax optimization, how to verify the global opti- 
mality of its solution, and the simplification utilizing the 



Here % is a Hilbert space that includes the range of 
p, and any Hermitian operator X satisfies the non- 
overestimating condition of Tr(Xpi) < £{p\) Vpi € % 
so that Tr(Xp) is a strict lower bound of £{p). X be- 
haves as a generalized entanglement witness operator: In 
the convex structure of quantum states, the hyperplane 
of po satisfying Tr(Xp ) = £ isolates the convex set of 
{pi | £(pi) < £q}, namely, Tr(Xpi) < £q for any p\ with 
£{pi) < £o- Hence X "witnesses" entanglement of finite 
amount £ , as Tr(Xp) > £ Q implies that £{p) > £q. Note 
that the conventional entanglement witness P3] detects 
entanglement as it isolates the set of separable states with 
£q = 0; for mathematical simplicity, in Eq. (2]), we use 
the different convention of a witness operator by the sign 
factor of —1 from the definition in literature. 

The optimal entanglement witness operator [TDHH1IH] 
X p satisfies Tr(X p p) = £(p), quantifying entanglement. 
It depends on the target state p and the choice off. From 
the viewpoint of optimization theory, X — X p is the glob- 
ally optimal point in that maximizes Tr(Xp). The 
global optimum is not necessarily unique for a given p, 
and any single global optimum suffices to obtain £{p). In 
some cases, X p does not exist in M%, and the supremum 
in Eq. ^ is only asymptotically accessed by an asymp- 
totic form of a witness operator [16]. The existence of 
X p depends on the choice of %. 

We mention one important finding of Ref. [16) . The 
duality between Eqs. ([I]) and ^ affords the criterion 
that is used to verify the global optimality of a solution 
of Eq. p}. From the fact that J^iPii^i^Hi) < £ {p) < 
'Yl, i Pi£\X^i)) 1 where {pi, \tpi}} is a pure state decomposi- 
tion of p = X^-Pi 71 "^' we nave the following equivalence: 
If X is the optimal witness operator of p, p belongs to the 
convex hull of Fx, convP^. The converse is also true, 



X = X p p e conv F x 



(3) 



symmetries. In Sec. IV we apply our approach to exam- 
ples, and discuss some general properties of the form of 
the optimal witness operator. We conclude in Sec. [V] 



where Fx is the pure-state set of {IV') | Tr(X7r^) = 
£(\ip))y. This implies that we can verify whether a wit- 
ness operator X is the optimal one for p, by checking 
whether p is decomposed into a convex combination of 
the elements of Fx ■ Hence, it is possible to check whether 
a numerical output of the optimization of Eq. ^ is at a 
local optimum (giving only a lower bound of £ ) or at the 
global optimum (with the exact value of £). We provide 
a quantitative way of checking this in Sec. |III C| 



II. OPTIMAL WITNESS OPERATOR 



III. MINIMAX OPTIMIZATION OF WITNESS 



We sketch how to quantify mixed-state entanglement 
by the optimal witness operator, which was developed 
in Ref. [16] . The primal problem in Eq. ([T]) is equiva- 
lently rephrased as the dual problem [TUlll6] that £{p) is 
obtained by the optimal operator X p , 



£(p)=T I (X p p)= sup Tr(A», 



:{X|Tr(X7r^)<£(|V>>) V|-0) e %}. 



(2) 



This section contains the followings. In Sec. |III A 



formulate the minimax optimization of X p . In Sec. IIIB 



we devise new techniques to overcome the difficulties that 
arise in the minimax optimization. In Sec. |III C[ based 
upon the correspondence of Eq. ([3]), we give how to quan- 
titatively check whether the optimization result is glob- 
ally optimal. In Sec. IMP] we show that our optimization 
is simplified when £ and p have common symmetries. We 



state the algorithm for the optimization in Sec. Ill E 



3 



A. Formulation 

To obtain X p , one optimizes the parameters of X un- 
der the condition of Ti(Xn^,) < 8(\ip)) in Eq. ((jj). To 
handle this condition efficiently, we translate it into an- 
other optimization problem. This idea of "optimization 
of optimization" leads to the essential part of our theory: 

Theorem 1. The optimization of the optimal witness 
operator X p in Eq. ^ is formulated as 



8{p) = Tr(A» 



sup inf F(IL|V)) 

n >olV'>eW 



sup GOT), 

n>o 



F(U, = Tr((n - pl)p), 

/i(n > |^)) = r &(n7r^)-f(|v» 1 

G(U) = Ti((U- m l)p), Mn = 



sup /i(n, 

Wen 



(4) 



Proof. We decompose each Hermitian operator X e H 
into X = II — fil with II > 0. Since X does not over- 
estimate 8, p is lower-bounded by /i > /in- For a given 
II, the finest witness operator is X = II — pul since a 
smaller p gives a finer X. Then, the expectation value 
TrXp is optimized by varying II, leading to Eq. Q. □ 

The minimization ("inf") directly reflects the non- 
overestimating condition in Eq. We note that Xp 
has an asymptotic form, when the supremum in Eq. (H) 
is asymptotically accessed; in this access, Tr((II — finl)p) 
converges, while II and pn diverge. 

The optimization in Eq. Q is a kind of the linear 
minimax problem [17j . The minimax problem has been 
widely applied to n-person games, finance, economics, 
and policy optimization. Since it is an optimization prob- 
lem ("sup") that includes a sub-optimization ("inf"), it 
is non-trivial to solve due to the following difficulty. In 
a simple-minded approach, one minimizes F(H, \ip}) over 
\ip) at every fixed II, and then maximizes the resulting 
G(II) over II. To apply usual competitive optimization 
algorithms to the maximization of G(n), one needs to 
evaluate the first or the higher order derivatives of G(II). 
The evaluation of the derivatives requires auxiliary min- 
imization to attain G(II + dH) for small deviation dH in 
every direction in the domain of II. Then, the computa- 
tional cost increases enormously as the dimension of the 
parameter space grows. 

Fortunately, there exists an efficient strategy that re- 
duces the computational cost in evaluating the deriva- 
tives of G(II). We here sketch the strategy following Ref. 
[17] . First we define the solution set E(II) 



i(n) = {\ip) €n\F(n,\ip)) = G(iT)}, 



(5) 



which becomes equivalent to Vx when X = II — pul- 
Namely, the elements in K(II) provide the value of G at 
II. The set K plays a key role in obtaining the derivatives 
of G; it is also useful for checking the global optimality 
of the solution of Eq. Q as shown in Sec. IIIC The 



behavior of G near II is governed by R(n); for small 
deviation dtt, G(II + dll) = inf|^ )eR(n) F(Il + dll, \ip)). 

When R has multiple elements, the multivariate func- 
tion G is not differentiable, i.e. the derivative of G cannot 
be expressed by well-defined Jacobian or gradient. For 
this case, we use DsnG(H), the directional derivative of 
G along SH at II, and D 5 jjG(II), that along the direction 

of steepest ascent 511 of G(II), 



D sn G(IL) = inf Tr [(p - 7r^Trp)<m] 

|V>eK(n) 

D 5(l G(Il)= sup D sn G(n) 
(5ner(n) 

ll<5ni|=i 



(6) 



where r(II) is the cone of possible directions of SH at II, 
and the directions are normalized, \\5H\\ = \\5Il\\ = 1, by 
the norm of ||A/|| = VTtM^M. 

Utilizing Eq. @, G(IT) can be maximized by the steep- 
est descent method. However, there are still two difficul- 
ties (i) that the direct evaluation of the steepest ascent 
using Eq. Q is another minimax problem, which may 
be laborious, and (ii) that the Hessian (the second or- 
der derivative) of G(n) may be unavailable when G(n) is 
not differentiable (i.e. when M.(H) has multiple elements). 
Note that the Hessian is used in most of efficient opti- 
mization algorithms since it informs the local behavior of 
the function to be optimized, which cannot be obtained 
from the gradient. For example, the Hessian determines 
whether a locally optimal solution is stable or unstable 
by its eigenvalues. We below devise new techniques to 
overcome the above difficulties. 



B. Minimization of F(I1, \ip)) 

To solve the minimax problem in Eq. Q, we need 
to accurately solve its "inf" problem, avoiding the con- 
vergence to local optima. Below, we propose the tech- 
niques that improve the optimization performance in the 
inf problem. We emphasize that the techniques are ap- 
plicable to general minimax problems. 

Parallelization. — The inf problem can be solved by 
adapting a parallelization approach with multiple solvers. 
A solver means a single optimization procedure that finds 
an (local or global) optimum, following a certain algo- 
rithm and starting from an initial guess (randomized or 
not) of the solution. We collect the results of the solvers, 
and assign the best value among them to G. Then M(n) 
is comprised of |-0)'s providing the best value or the suf- 
ficiently good values slightly deviated from the best one 
within a given tolerance. M.(H) can be fully identified, in 
principle, given that the unlimited number of solvers are 
available. In practical situations, however, parallelization 
details such as the number of solvers and initial guesses 
need to be tuned depending on the problem. 

To design efficient parallelization, one may require to 
investigate entanglement classes with respect to a given 
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measure 8. The entanglement class is an exclusive sub- 
set of quantum states whose elements share the common 
type of entanglement. For example, there are two entan- 
glement classes for two qubits: one is entangled and the 
other is separable. When confined to pure states, some 
classes, e.g. two-qubit separable states, have vanishing 
volume. It is generally possible that some elements of 
K(n) belong to different entanglement classes, even to 
classes with zero volume. However, if a solver is designed 
to minimize ^(n, over the full set of pure states at 
once, the zero-volume classes can be rarely accessed and 
the solution tends to appear in classes with hnite volume. 
This can lead to the insufficient construction of M(II). 

This problem is resolved by the divide-and-conquer ap- 
proach: (i) divide the full set of pure states into the en- 
tanglement classes with respect to the target entangle- 
ment £ , (ii) assign solvers to different classes, and then 
(iii) let each individual solver minimize F(H, over 
\tp) only within the assigned class. In this approach, all 
the classes can be thoroughly accessed. Though the full 
classification of entanglement for more than three qubits 
is complicated, even an incomplete classification, such as 
the dichotomy between finite and vanishing £, improves 
the performance of the optimization. 

Finding elements o/R(n). — We need sufficient infor- 
mation of R, to precisely evaluate the directional deriva- 
tive of G or to verify the global optimality of the opti- 
mization result. R may be numerically found by consid- 
ering the local optima of the inf problem, as mentioned 
above. This procedure has two non-trivial points: (i) 
The values of F at the local optima may deviate due to 
numerical errors, and (ii) all local optima residing in the 
high-dimensional domain are hard to be found. 

We handle the point (i) by extending the concept of 
the (exact) solution set R to the set R of the candidates 
for the solutions. R consists of not only the global op- 
tima but also the local optima which give the value of 
F sufficiently close to the globally optimal one. This ex- 
tension makes G more robust against errors in numerical 
computation, avoiding the possibility that the global op- 
tima are treated, due to the errors, as local optima and 
excluded from R; the exclusion will induce the failure in 
evaluating D$jjG, hence the maximization of G. 

On the point (ii) , the parallelization of the inf problem 
is helpful, as it provides a number of local optima, the 
elements of R. Sometimes, more elements of R are re- 
quired. They are provided by the local maximum points 
\ip) E H of the function 



d(ivo jv>o)) - Mm iv>>) - m iv>o>)) 



(7) 



where d(\4>), \ipo)) measures the distance between l^) and 
the best optima IV'o) obtained from the parallelization, 
and k is a positive number fixed arbitrarily large. The 
exploration of the local maximum points by Eq. Q de- 
picts the search of \tp)'s that are located far from the best 
known optimum l^o) but provide the value of F close to 
the best one; the second term in the cost function of 
Eq. Q acts as the barrier that rules out the points 



providing not-small-enough values of F. With reason- 
able tolerance for F values, the above search accesses all 
available candidates of R. The search can be generalized 
by considering also the distance from the other locally 
optimal candidates in R. Note that if a better optimum 
point giving a smaller value of F than |^o) appears in 
this search, then |^ ) should be replaced with this bet- 
ter one. This search complements the parallelization in 
finding the global optima of the inf problem. 

Smoothing of G(n). — The function G resulting from 
the inf problem can be non-analytic, as mentioned in Sec. 
Ill A To surmount the non-analyticity and to efficiently 
evaluate the derivatives of G, we approximate G to an 
analytic function G within good accuracy. We call this 
approximation the "smoothing" of G. 

In this method, we use the candidate set R = 
{1-01), • • • , IV'Af)} instead of R, as the computation of G 
becomes more robust against numerical errors. When R 
has enough elements, G is faithfully approximated by G, 

GOT) = H(F(U, \rp 1 )), • • • , F(n, |^»; b). (8) 

H(oex, a%, • • • ,aN]b) is an analytic function imitating 
min{a„} with a smoothing parameter b > 0, satisfying 

H(ai, • • • , oln] b) < lim H(ai, • • • , oln] b) = min a n . 

b— >0 n=l,— ,N 

Larger b makes H smoother, but more deviated from 
min n= i t ... ^{a n }. As G(II) is a differentiable function 
with gradient VG, the directional derivative in Eq. Q 
is approximately obtained by a simplified form 



AmG(ri) = Tr(VGflT 

dH(n,{\^)}) 



(9) 



VG = 



an 



Similarly, the direction of the steepest ascent is given by 
VG/|)VG||. Accordingly, we employ the gradient and 
the Hessian to maximize G efficiently and precisely, with 
avoiding the difficulty in directly evaluating Eq. (|6|. 
We mention our choice of the smoothing function H, 



H(ai, ■ ■ ■ ,a N ; b) = H(H(a 1 ,- ■ ■ , a N -i; b), a N ; b). 

For N > 3, H is not exactly symmetric under the per- 
mutation of a n , but, the asymmetry is negligible when b 
is small enough. Then Eqs. ^ and ^ arc rewritten as 

G(n) = Tr (lip) - /in, 

An = -H(-n(U, |V>i}), • • • , -M n > l^iv)); b) > Mn, 
dfi u 



VG = p 



an ' 



Here, we have used the property of H{a\ + a', • • ■ , ajv 
a'; b) = H(a\, • • • , an; b) + a'. 
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C. Verification of global optimality 

The maximization of G in Eq. Q provides a candidate 
of the optimal witness operator. We will show that in our 
approach, one can verify whether the candidate is the 
global optimum of the problem, i.e. the correct optimal 
witness operator. 

We first translate the duality correspondence in Eq. 
|3]) into the following, to use it in the minimax approach. 

Theorem 2. A witness operator X = II— finl is optimal 
for £(p) if and only if d m i a (p; convR(n)) = 0. Here, 
rfmin(p';S) = min CT gs|jp' — ct||hs measures the minimal 
distance between a state p' and a set S of states a by the 
Hilbert- Schmidt norm || • • • ||hs- 

Proof. Px is equivalent to R(II) when X = IT — /inP And 
p G convPx if and only if d m i n (p; convPx) =0. □ 

Based on Theorem [2j we use d,„i n (p; convIR(n)) for 
quantitatively judging how far a candidate II for the so- 
lution of Eq. Q is deviated from the global optimum. 
Concretely, d m i n (p] convR(IT)) is not a distance between 
the candidate and the global optimum, but a degree of 
the deviation from the duality correspondence in Eq. ([3]) 
which should be satisfied by the global optimum. Hence 
any a priori information about the global optimum is not 
necessary, and the only requirement is that the elements 
of R(n) are sufficiently found by the parallelization of 
the inf problem or by the exploration using Eq. 0. 

To use Theorem [2] in numerical optimization, the ex- 
tension of 1 to 1 and that of d m i n (p; convR(IT)) to 
d m i n (p; conv R(n)) are crucial. When some elements of R 
(therefore those of conv R (II)) are accidentally excluded 
due to numerical error, d m { n (p; convR(IT)) can be over- 
estimated even if II is the globally optimal one. On 
the other hand, given an appropriate value of the toler- 
ance for F(n, - G(IT), R(II) is the set of all relevant 
|V>)'s. The tolerance judges whether a state \ip) is rel- 
evant to M(II) by considering how accurately £ (\ip)) is 
quantified by X = U — p,nl, since £(\tp}) — Tr(X7r^) = 
F(II, - G(IT) because of Trw^ = Trp = 1. 

The choice of the tolerance for F(IL, \ip)) — G(H) gov- 
erns the credibility of d m i n (p; convR(IT)). Too small tol- 
erance is not helpful, as it causes the problem similar 
to the case of R. Too large tolerance causes the inclu- 
sion of irrelevant elements |V')'s in convR(n), for which 
Tv^Xtt^) is far deviated from £(\ip)). It leads to the 
underestimation of d m i n (p; conv R (IT)), resulting in mis- 
leading verific ation of the global optimality. In our ex- 
amples in Sec. 



IV 



the values of F(n, \ip}) at \ip) G 
are well-separated from the values at other local optima, 
hence the adequate value of the tolerance is straightfor- 
wardly chosen. 

The result of d mm (p; conv R (IT)) = docs not only 
confirm the global optimality of witness X = II — pul, 
but also provide the optimal pure-state decomposition 
of p for £; a convex sum € convR(II) is the optimal 
decomposition of p, when it is identical to p and has the 



vanishing minimal distance of d m i n = 0. It is because 
any convex sum p' G convR(II) = convP^ itself is the 
optimal pure-state decomposition by Theorem 2 in Ref. 

We note the possibility that at some IPs, G(IT)'s are 
similar to each other, but d m i n (p; convR(IT))'s are largely 
different; see the example in Fig. [2] This indicates that 
d m in{p', convR(II)) is a more useful tool for checking the 
global optimality than the mere comparison of G(IT)'s. 



D. Simplification by symmetry 

The common symmetries of £ and p simplify the min- 
imax approach. It is identified by the symmetry group 
G = {U I W = If- 1 , p = UpU\ £{p') = £{Up'W) Vp' e 
H}. Usually G has local unitary operations and permu- 
tations of subsystems; see Sec. |IV| 

We first introduce some notions for later use. We call 
O G H a symmetric operator with respect to G when 
O = UOW holds VU e G. We also call C = {P l eH} & 
"basis" set of symmetric operators Pi with respect to G, 
when C consists of the maximal number of linearly in- 



dependent P^s satisfying Tr(PiPj) 



Any symmet- 



ric operator O is decomposed as a linear combination of 
Pi's, and the coefficients of the decomposition can serve 
as the parametrization of O. And, any operator O G % 
is symmetrized as O s = £ 4 Tr(OP) P = 'EueG UOW; 
it is obvious that O = O if and only if O is symmetric. 
Note that though the choice of C may not be unique, any 
choice provides the same result of the symmetrization. 

In addition, we call a set S a symmetric set with respect 
to G when § = {UOW \ VO G §, Vf7 G G}; accordingly, 
a set of symmetric operators is symmetric. There exists 
a minimal subset S A of a symmetric set S such that S 
is generated from S A by applying the elements of G to 
the elements of § A . We call S A the asymmetric unit of 
S. It contains the information of S not related to the 
symmetries. The choice of S A may not be unique, but 
any choice generates the same S. 

Based upon the notions, we find the following features: 

Lemma 1. There exist the optimal witness operator X p 
and the corresponding Tl p symmetric with respect to G. 

Proof. Consider the operator of X = X p U for an op- 
timal witness operator X p and any U G G. Then A" is a 
witness operator X G M% not overestimating £ since 
TriWXpU-K^ = Tr(A p £/7r [/t) < £(U\i/>)) = £{\^)) 
holds VJ7 G G and V|^) G %. The witness X is optimal 
for p, since £(p) = Ti(X p p) = Tr(Ap). Hence the sym- 
metrization X p = Tr(ApPi) Pi of a non-symmetric 
optimal witness operator X p constructs the optimal wit- 
ness operator X^ symmetric under VU G G. When X p is 
symmetric, so is n p because I is symmetric. □ 

Lemma 2. If X and n are symmetric with respect to G, 
Fx and R(n) are also symmetric. 
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Proof. V|V> G fx and VU G G, £(\i/j)) = Tr(XTT^) = 
£{U\if))) = Tr(WXUir^). This shows that J7|^) G P x , 
indicating that P^ is symmetric. And R(II) is equivalent 
to P x for X = n - fj, n I. □ 

From these two Lemmas, we derive the simplification 
of the minimax optimization of the optimal witness op- 
erator, which is the main result of this subsection: 

Theorem 3. Suppose £ and p have the common sym- 
metry given by G. Then in solving Eq. the domain 
of the sup problem II > and the solution set of the inf 
problem M(II) can be replaced by II s > and [M(II S )] A 
( with respect to G ), respectively. 

Proof. We prove this theorem by the following two parts. 

(i) n -> n s : Lemma [j] ensures the existence of n o 
symmetric with respect to G. Hence the optimum £{p) 
is accessible in the sub-domain of symmetric operators 
II s > 0. It proves the replacement of II by II s . 

(ii) R(n) -> [M(n s )] A : We start with M(li s ) that is sym- 
metric with respect to G, replacing II by II s as in the part 



(i). We consider a pure state 



'). With small 



enough r > 0, one defines B r (\ip}), an open ball of radius 
r centered at here the set of pure states is treated 
as a metric space with a Hilbert-Schmidt norm metric. 
Then for any U G G, the set of {U\b) | V|6) G B r (\ip))} be- 
comes B r (U\ip}). Therefore! = B r (|^i))UB r (|V) 2 })U- • • 
for IR(II S ) = 1^2): • • • } is symmetric with respect 

to G. Due to the symmetry, {P(II S , \ip)) \ G B} = 
{P(LI S , \ip)) I € B A } holds. It implies that the behav- 
ior of G near n s is fully described only by [R(II S )] A . □ 

The replacement in Theorem [3] reduces the computation 
cost, since a symmetric operator is parametrized by a 
smaller number of parameters than a non-symmetric one. 
Namely, the dimension of the domains for the optimiza- 
tion is reduced by the symmetry. This reduction makes 
the optimization easier to converge to the global optima, 
since local optima (which are deep enough to be stuck) 
tend to appear less often in a lower dimensional space. 

The minimal distance in Theorem [2j the quantitative 
criterion to verify the global optimality of a witness oper- 
ator, is also streamlined by exploiting the symmetry, i.e. 
by replacing d m i n (p; convR(II)) to d m i n (p; conv]R(II s )). 

We below apply the symmetries to numerical proce- 
dures. Using a "basis" set C = {Pi}, we map operators 
in H to real vectors in Euclidean space. Substituting 

n s = J2i v i p * and p = Ei r i p i into E q- ©' we § et 

£(p) = sup inf (v,r}-«v,q}-£(|V)))Trp, 
n>o 

where (•, •) is an inner product of vectors, [v]j = Vi, [r]; = 
r h and [q\i = Tr(P l7 r ). Rewriting F(U, \j>)) = P(v,q) 
and G(II) = G(v), the directional derivative of G at v in 
direction v' (||v'|| = 1) becomes [see Eq. ([6])] 

D V ,G= inf (r-Trpq,v'). 

|v>)e[R(ns)]A 



To smooth G, we find the elements of [M(II S )] A as the 
local maximum points \ip) € Ti. of the function in Eq. 

0. 

llq-qoll+MM^MV))-^ 8 ,!^})), 

where p,(U s , \ip Q )) = p, n s, [q ]j = Tr(P;7r^ ), and k 
is a large positive constant. Then the gradient of the 
smoothed G [see Eq. &\ is 



VG 



dH(-(i(v,qi 



-/i(v,qjv);6) 



<9v 



where p(y, q*) = p(Il, \ipi)) = (v, q.i) -£(\ipi)) for [q j; ] 3 = 
Tk(PfTr^) and [i(II s )] A = {|^i), • ■ ■ , |^»- 

The minimal distance is also rewritten from the 
Hilbert-Schmidt norm to the Euclidean norm of 
rfmin(r';S) = min se s ||r' — s||. Then the sufficient and 
necessary condition of the global optimality becomes 



(2min(r;conv{qi,--- ,qAr}) = 0, 



(10) 



where [r], = Tr (Pi p), [q^ = Tr(Pj7r^, 4 ), and 
{\tpi), • • • , \?Pn)} = [M(II S )] A for a non-smoothed G and 
• • • = [§(n s )] A for the smoothed G. 

By finding r' = £). w-q* G conv {qi, • • • , q^} that pro- 
vides d m - m = 0, we build the optimal pure-state decom- 
position p=J2z w i^: where n % = Ec/eG Un^W. We 
note that [M(n s )] A or [§(n s )] A holds the key informa- 
tion of the convex structure of ]R(II S ) with respect to 
£ , hence, is the minimal set for analysis, since K(II S ) is 
constructed from [M(II S )] A by the symmetry operations. 



E. Algorithm 

We summarize this section outlining the entire algo- 
rithm for the minimax optimization of witness operator 
for a general state p and a measure £ . The steps of the 
algorithm are as follows. 



1 



Investigate the common symmetry of £ and p, and 
construct the corresponding "basis" set C = {Pi}. 
This step reduces the number of independent pa- 
rameters. [Sec. 



HID 



2. Inspect the entanglement classes of pure states re- 
garding £ and divide the set of pure states into the 
classes. It improves the convergence of the prob- 
[Sec. ~ 



lem. 



IIIB 



3. Set the number of parallel solvers and define the 



smoothing function. [Sec. Ill B 



4. Solve the minimax optimization; maximize G(II) 
using a competitive optimization algorithm. At 
each II, the value and the derivative of G(II) are 
obtained via the parallelization and the smoothing. 



[Sec. IIIB 
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5. Verify the global optimality by evaluating 
d m \ n {p\ convR(II)). [Sec. Ill C| When the solution 
is far from the global optimum, increase the 
number of parallel solvers and/or subdivide the 
entanglement classes. Then retry the optimization. 



IV. THREE-QUBIT GHZ ENTANGLEMENT 

We apply the minimax optimization to the quantifi- 
cation of the three-qubit GHZ entanglement in example 
states. The example states are composed of three species: 
The GHZ state |GHZ) = ^ (|000) + | 111)), the W state 

\W) = ^(|100) + |010) + |001)), and the white noise 1/8. 
We reproduce the result of Ref. [TB] for the mixtures of 
two out of the three species, and scrutinize the mixture 
of the three species. The examples reveal some general 
properties of the form of the optimal witness operator 
and of the optimal pure-state decomposition, hence, of 
the convex structure of general mixed states. 

This section consists of the introduction of a mea- 
sure of three-qubit GHZ mixed-state entanglement in 
Sec. |IV A[ how to apply our minimax approach to the 
measure in Sec. |IVB| three examples of quantification of 



IV E and 



GHZ mixed-state entanglement in Sees. IV C 
the usefulness of <i m i n in the numerical quantification in 
Scc. |IVF[ It also has Sec. |IV G"] where some general prop- 
erties of the form of the optimal witness operator are 
discussed. 



A. Extensive three-tangle 

The entanglement structure of three qubits is more 
complicated than that of two qubits, as there are six 
classes in the former while only two classes (entangled or 
not) in the latter. Among the six classes, there are two 
genuine tripartite entanglements, GHZ and W, whose 
representative states are |GHZ) and \W) respectively 
[15] . The GHZ\W class, the set of states with finite 
GHZ entanglement, comprises the outermost part of the 
convex structure of three-qubit mixed states [2D] . 

We choose the extensive three-tangle T3 PS] as the 
measure of three-qubit GHZ entanglement. It quantifies 
the three-qubit GHZ entanglement and generalizes the 
three-tangle T3 [2T] to be invariant under stochastic local 
operations and classical communications (SLOCC) for 
mixed states [ini[22]. It is defined as "Tad ip)) = ^/r 3 (|V>)) 
for pure states, and extended for mixed states by the 
convex roof construction. The name "extensive" is given 
by the property T3 (rp) = rT 3 (p) for positive real r. The 
optimal witness operator for r 3 can be also constructed 
by Eq. Q, but we do not discuss it in the present work. 



B. Minimax optimization for 75 

When one applies the minimax formalism in Eq. Q to 
75, the domain of pure states for the inf problem needs to 
be carefully considered due to the W class. For a mixture 
p of GHZ- and W-class states, Fx p or R should include 
the states in both the classes. However, the set of W-class 
states has zero volume with respect to pure states, which 
means that the W-class states are hard to be accessed 
when a pure state is parametrized to cover the overall 
pure states. Even more, the directional derivative of 
73(1^)) at the W-class pure states in the direction along 
the GHZ-class pure states shows divergence, originated 
from the definition of T3; D\^% = (2^/t^) -1 D^t 3 di- 
verges when 7jj = T3 = and D^t 3 is finite. This diver- 
gence obstructs the inf problem to be solved correctly. 

To avoid the above problems, we separate the domain 
of the inf problem into the GHZ\W and W classes, as 

£(p) = su P inf Tr((II-£i», 

n>o |^G>eGHZ\w 

\lP w )£W 

p = H(p w , p G ;b) > p = max{p G ,p w }, 
p G = Tr(n7r^ G ) - £(|V>g)), 
p w = Tr(n7r^ w ). 

Here we parametrize \if)w) an d \ipc) by using the gener- 
alized Schmidt decomposition [23], 



|V) = Ui <g> U 2 <8 U a \V), 
= Ao|000) + Aie**|100) 
+ A4IIII), 



A 2 |101) +A 3 |110) 



where f/j's are local unitary, A; > 0, and Af = 1- 
T 3 (\4>}) = 75(|*» = 2A A 4 holds for GHZ\W-class states 
and A4 = 4> — for W-class states. 

We mention the optimization settings chosen for the 
examples below. We parallelize the inf problem with 30 
solvers. Every solver in the sup and inf problems runs 
based on the sequential quadratic programming [M] al- 
gorithm, an iterative method for nonlinear optimization. 



C. Mixture of |GHZ) and 7/8 

First our formal approach is applied to the mixed state 
of |GHZ) and the white noise 1/8, pgi(i) = (1 — q)^GHZ + 
qI/8 with q £ (0,1). The entanglement quantification 
for this state was analytically studied in Ref. [TB]. We 
reproduce the result of Ref. [TB] numerically. 

The symmetry group for pci consists of (i) the local 
phase rotation 

R h (a,p) = (7r + e M ^i)(»K + e l,9 7r 1 )®(7ro + e-' i ( a+ ' 3 )7r 1 ) 

for the continuous range a,j3 € [0, 2ir), (ii) all permu- 
tations among three parties, and (iii) 0-1 flip for three 
parties erf 3 . Here tt — |0)(0|, 7Ti = |1)(1| and a x is a 



Pauli matrix. We select the bases of symmetric operator 
as Pi=oA = £j/v2 and P 2 = (I ~ Ho)/y/6, where we 
define £,=0,1.2,3 as 



£0 = "000 + Till, 

S 2 = -i|Q00)(lll| - 



Ei = |000><lllj - 
-»|111><0001, s 3 



|111>(000|, 
""000 — TTui, 



with ttooo = |000)(000| and 7r m = |111)(111|. 

We obtain %,{p G i(q< <7o)) = 1 - <?/<7o and %,{pgi{<1 > 
go)) = for g — 0.304, where the global optimality 



of X PGI is ensured by d m ; n (p; conv]R(n)) 
optimal witness is also found as 



10- 



The 



X 



PGI 



4.053P + l-604Pj - 3.000/ for q < q , (11) 



and X PG1 — for q > 0. The optimal decomposition for 
q < q is p G i{q) = 7i(g)7r GH z + (1 - %{q))pz, where p z 
is the W-class mixed state 



Pz 



1 

18 



E 



s=± 
ro,n=Q,l,2 



(0 e l2 T ) 



|Z s oo)i 



l z +oo) = 2J a »+J"+fcl^' fc )' 
i,j,fe=0,l 

|Z_00) = 2J a 3-(i+i+fe)|y^) 5 
jj,fc=0,l 

and (ao, 01,02,03) ^ (0.7436,-0.1750,-0.2133,0.4677). 
For the above decomposition, the asymmetric unit of 
Fx has two elements. We choose them as |GHZ) and 
|Z+oo)j as the other |Z smn )'s are obtained from |Z + oo) by 
the elements in the symmetry group. 



D. Mixture of |GHZ) and \W) 

Next we study pGwip) — (1 — p)""ghz + P^w with 
p S (0, 1), and numerically confirm the analytic result in 
Ref. |16) . The form of the optimal witness operator for 
Pgw depends on the Hilbcrt space H under considera- 
tion. It has an exact form in the range of pgw, while 
it has an asymptotic form in the full Hilbcrt space of 
three qubits. Both the forms give the same result of 
T 3 (p G w(p < Po)) = 1 ~p/Po, %(pgw(p > Po)) = for 
Po ~ 0.3731, and the optimal pure-state decomposition 

Pgw(p) = 

/3(P)7TGHZ + 3 2^ ""Z'nCfo)' P<P0, 

5 E ^; 



n=0,l,2 



P>P0, 



where \Z' n (p)) = yT^|GHZ) - 
class states with Tr(X PGW 7rz^) 



^tp e 2"«/3|w) are W- 
= 0. Without loss of 
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FIG. 1. (Color online) Extensive three-tangle T3 of the 
three-qubit full-rank mixed state of pGwi{p,q) = (1 — P — 
<?)ttghz +pKw + q!/9>- The physical region of (p, q) is defined 
by p, q, (1 — p— q) g [0, 1]. The states on the abscissa and on 
the ordinate are pgw(p) and PGi(g) whose values of 75 vanish 
for p > po — 0.373 and q > qo ~ 0.304, respectively. 



generality, Px PGW ] consists of |GHZ) and |Z (po)) in 
the case of p < po, while a single element \Z' (p)} for 
p > Po; the others |Z^(p))'s with n = 1,2 are obtained 
from |Z' (p )) by the symmetry operator of Pl(^P, 

By using this example of Pgw, we show how the 
asymptotic form of the optimal witness operator is han- 
dled in the minimax formulation. In this case where H is 
the full Hilbcrt space of three qubits, the elements of G 
are (i) discrete local phase rotations Pl(^t, ^n) with 
n = 0, 1, 2, and (ii) all possible permutations among the 
three parties. Note that the symmetry group for Pgw 
is the subgroup of that for pgi- The bases of symmetric 
operator are 



-Pi=0,l,2,3 — Si/V2, 
Pi = TTWi , P5 = {^W 2 + n W 3 ) /2, 

A 



(12) 



Pe = tt Wi , P 7 = (tt 



w 2 



'w 3 



where \W n ) = R h (0, 2 -^^)\W) and \W n ) = a r f\W n ). 
To tame the divergent parameters of the optimal wit- 
ness operator of an asymptotic form, we set the bound 
of the parameters of n = '}2 li v i P i as \vi\ < k where 
k is an arbitrarily large positive number. The accu- 
racy of the result increases as k grows; for k = 10 2 , 
both of d min (p; convM(n)) and 73(pgw)-^P^gwPgw) 
are ~ 10~ 2 , and for k = 10 3 they both decrease to 



10" 



We obtain, for p < p , that Tr(X PGW Pj = 



=0,1 



1/V2, Trpf PGW P 4 ) -> 1 - po" 1 , and the other 
Tr(X PGW Pj^o,i,4)'s diverge, as k increases. Regardless 
of the divergence, Tr(pGW-^"p G w) successfully converges 
to the exact value of £ (pgw) since Tr(pGwP^o,i.4) = 0. 
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E. Mixture of |GHZ), \W), and J/8 

We consider a more complicated example of the mix- 
ture of | GHZ), \W) and 1/8, p GW i{p,q) = (1 - P - 
9)ttghz +qI/& with p,q,p+q G (0,1). The symme- 

try of pgwi is the same as that of pew! see Eq. jl2|). In 



Fig. [I] we numerically compute 73(pgwi) using 
imax optimization. The global optimality of the optimal 
witness operator is ensured by d m i n (p; convIR(n)) ~ 10~ 4 
in the computation. 

We address the optimal pure-state decomposition 
of pgwi (Pj q) f° r the case that pgwi (p, ?) has finite 
7a(PGWi(P) 3))- The decomposition form is numerically 
obtained, after finding the optimal witness operator. It 
consists of one GHZ\W-class pure state and multiple W- 
class pure states as pgwi = a7T r + (1 — ol)pz, where 
a = a(p,q), \r) = vT^IOOO) + y/r\lll) is a GHZ\W- 
class state, r = r(p,q) € (0, 1), and pz is a mixture of 
W-class states; this decomposition form indic ates that 
% has the form of 71s(pgwi) = OL%(\r)) = 2a^/r{\ - r). 
pz is decomposed into pure states in the symmetric set 
whose asymmetric unit has two elements \x) and 

Pz = PlPx + p2Py, A+/3 2 = l, Pi,fo>0, (13) 

1 2 

PV = \ E *Vm», ^^RLi 2 -^, 2 -^)^), 
n=0,l,2 

i,j,k=0,l 

\y+) = c o |000) - Cl \W 2 ) - c 2 |W 3 ) + c 3 |lll), 
|y_> = c |000) - ci|W 3 ) - c 2 |W 2 ) + c 3 |lll), 

where bj's, Cj's and /9j's depend on p and q. 

We describe how the optimal decomposition of pgwi 
is related to those of pci and pew- For pGWi with finite 
75, [Px p ] A consists of one GHZ\ W-class state \r) and 
two W-class states \x) and On the other hand, 

[Px p ] A contains one GHZ\ W-class state |GHZ) and one 
W-class state |Z +0 o) for pci with nonzero % (i.e. q < 
qo), while one GHZ\W-class state |GHZ) and another 
W-class state |Zg(p)) for p GW with nonzero % (i.e. p < 
po). As pcwi(p,q) continuously approaches to pgi(<z) 
by p — > 0, we see that \y + ) — > |Z +0 o), Pi — > 0, and 



02 



1. And, 



|Z'q>, A 



1, and /3 2 — > 0, as 



PGwxip, 9^0)^ Pgw(p)- 

This example of pcwi(p,q) is the simplest full-rank 
mixed-state that is composed of pure states with different 
genuine multipartite entanglement such as GHZ and W. 
Entanglement in this type of mixed states has not been 
quantified before, hence, this example demonstrates the 
usefulness of our numerical approach. As discussed in 
the following subsections, this example demonstrates the 
usefulness of d m ; n in Sec. |IVF| and reveal some general 
properties of the form of the optimal witness operator 
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FIG. 2. (Color online) Quantification of entanglement in 
p(p) = PGWi{p,q = 0.0380) by 75. The exact value of Ti{p), 
Tr(Xpp), is numerically obtained by the optimal witness op- 
erator X p , and plotted by green filled circles (left ordinate) as 
a function of p € [0.01, 0.05]. The global optimality of X p is 
verified by the minimal distance d m in(p; convR(IT)) ~ 10 -4 , 
which is also plotted by purple open circles (right ordinate). 
This exact quantification is compared with the underestima- 
tion by a misleading less optimal (locally optimal) witness 
operator Xi cs . The underestimated value of 75 by Tr(Xi es p) 
is plotted by blue filled rectangles. We also plot the minimal 
distance d m i n for this underestimation case by red open rect- 
angles, which clearly shows that X\ cs is not globally optimal. 



and of the optimal pure-state decomposition, hence, of 
the convex structure of mixed states in Sec. IIVGI 



Usefulness of d m i n 



The example in Sec. |IVE| shows that the minimal 
distance d m i n (p; conv]R(n)) is essential for verifying the 
global optimality. In Fig. [2] the exact computation of 
T3 by the optimal witness operator X p (p) and a slight 
underestimation by a misleading locally optimal witness 
operator X\ cs {p) are compared for PGWi(p,q = 0.0380); 
the parameters of the two operators are compared in Fig. 
[3] The global optimality of the two operators are tested 
by the minimal distances d m i n (p; convM(n)). The values 
of Ti(Xp) by the two witnesses X p and Xi es differ by a 
small amount (~ 10~ 4 ) near p = 0.01, hence, it is diffi- 
cult to check the global optimality of a witness operator 
solely by the expectation vale of Tr(Ap). In contrast, 
<i m i n is clearly distinguishable by ~ 0.01 between the two 
witnesses, even in the case that the values of Tr(Ap) are 
very similar. This demonstrates the usefulness of d m ; n 
for distinguishing the globally optimal witness operator 
from other locally optimal (less optimal) ones, namely 
for testing the global optimality of a witness operator. 

We discuss the appearance of a misleading less optimal 
witness. When one searches for X p for general p and £ , 
a misleading less optimal witness A m may frequently ap- 
pear if the target state p neighbors a state Ph with high 
symmetry. In this case, X m shares the high symmetry 
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FIG. 3. (Color online) Parameters of II — Yli Vi ^ >i of the 
two witness operators in Fig. [2] the optimal witness operator 
X p (green filled circles) and the misleading locally optimal 
one Xies (blue open rectangles), which quantify GHZ entan- 
glement in Pgwi(p, 1 = 0.0380) by 75- We select the bases of 
symmetric components in Eq. (121, since pgwi has the same 
symmetry as pgw- 



with ph- Since the high symmetry results in many sym- 
metry operators, Px m has lots of elements generated by 
these operators, and accordingly convPx m exhibits large 
dimension and size. This indicates that X m is a local op- 
timum point to which the optimization easily converges; 
any small deviation of a witness operator from X m results 
in the decrease of Tr(Xp) since the deviation degrades the 
symmetry of witness and reduces conv Px m significantly. 
In the example of Pgwi m Sec. |IVE| the misleading wit- 
ness operator X\ cs has the form similar to X PG1 ; one can 
notice this by comparing the parameters of X\ cs shown in 
Fig. [3] with Eq. (11). In the parametrization space with 
respect to Eq. (12), the dimension of convPx POI is six 
while conv Fx 

^ D 



is only three-dimensional. Due to this 
difference, Xi es may appear frequently though it is mere 
local optimum. Such misleading local optimal witnesses 
can be ruled out by computing c? m i n - 



this approach the perturbative construction of the optimal 
witness operator. This idea is based upon the continuity 
of entanglement measure, \£{p') — £{p)\ — > as 6 — » 0, 
and upon the assumption that the form of the optimal 
witness operator or the form of the optimal pure-state 
decomposition changes continuously and smoothly from 
p to p' . This perturbative approach works well for the 
states inside each family of pci, Pgw, an d Pgwi- The 
success of the perturbative construction saves the number 
of parallelized solvers for the inf problem in Eq. Q . 

However, the perturbative construction is not always 
successful. For example, one may try to perturbatively 
construct the optimal witness operator for pcwi, start- 
ing from that for pgi or pewj i.e. , from the simpler 
case for which the optimal witness operator is known. 
It fails, because of the change of the symmetry (in the 
case from pgi) or the change of the rank (from pew)- In 
the case from pci, the witness for pgwi perturbatively 



constructed from X a 



As shown in Figs. p2] and 



[3] (see rectangles), Xi os (p) underestimates 73(pgwi(K <Z — 
0.0380)), is not globally optimal, and has a form differ- 
ent from X p (p), indicating the failure of the perturbative 
approach. In the other case from pgw> the parameters of 
the full-rank witness operator for pgwi diverge at pew 
(i.e. at q = 0); note that as discussed before, the optimal 
witness operator for pew has an asymptotic form, when 
one choose H as the full Hilbcrt space of three qubits. 

The failure of the perturbative construction is mani- 
fested by the non-smooth evolution of the convex struc- 
ture (with respect to £ ) from p to p'. The non-smooth 
evolution means that the optimal witness X p or the con- 
vex set convPx p evolves non-smoothly between the two 
cases. We notice this feature, by observing the form 
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of the optimal pure-state decomposition in Sees. IV C 
that Px p suddenly changes from pgwi to pci or 
n the case of pgi, Pjc p consists of infinitely many 
One of them is |GHZ) and the others are the 
l-Z+oo) by the symmetric opera- 



Pgw- 
states: 

states generated from 
tors of pgi, including the continuous local phase rotation 
Ri,(a,f3) with Va,(3 £ [0,27r). In contrast, there are ten 
elements £ Vx in the case of PGWii and only four for 
Pgw- Hence, the convex structure of convPx p is non- 
smoothly connected among pcwi; Pgi 7 an d pgw- 

This example shows that there is no single universal 
form of the optimal witness operator in the three-qubit 
case, in contrast to the two-qubit case [B], and that the 
convex set convPj5f p (or its evolution) needs to be studied 
to construct the optimal witness operator and to under- 
stand multipartite or high-dimensional entanglement. 



G. Perturbative construction of optimal witness 

When one constructs the optimal witness operator for 
p' = (1— S)p+Sa (S -C 1), one may guess that the optimal 
witness operator for p is useful, and try to construct the 
optimal witness operator for p' from that for p; we call 



V. CONCLUSION 

We develop a numerical optimization approach for 
finding the optimal witness operator that quantifies 
mixed-state entanglement. It is applicable to general en- 
tanglement measures and states. Importantly, the global 
optimality of the optimization result is verifiable. Further 
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development of our approach will provide a powerful tool such as nonlocality and quantum channels, 
for studying not only multipartite mixed-state entangle- 
ment but also other convexity-related physical properties 
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